#include <stdio.h>
#include <mpi.h>
#ifdef HAVEFITSIO
#include <fitsio.h>
#endif

#include "profile.h"
#include "kdtree.h"

#ifndef _HALOFINDER_
#define _HALOFINDER_

#ifdef MEMWATCH
#include "memwatch.h"
#endif

/* Some debugging macros
   undef DEBUG for no debugging
   DEBUG_LEVEL = 0 is for basic debugging
   DEBUG_LEVEL = 1 is for messages printed by a single task but not as critical
   DEBUG_LEVEL = 2 or above is used for messages that every task will print 
   
   define DEBUG_IO some output as follows
*/

#ifdef NDEBUG
#undef DEBUG
#define DEBUG_LEVEL -1
#undef DEBUG_IO
#endif

#ifdef DEBUG
#ifndef DEBUG_LEVEL
#define DEBUG_LEVEL 0
#endif
#endif

#define MAX_FILENAME 1024

#define BIG_G 4.304e-9
#define RHO_CRIT 2.77519737e11 /* Critial mass density  in h^2 M_sun/Mpc^3 */
#define CSOL 299792.458 /* velocity of light in km/s */

#define TAG_NUMDATA_BUFFEXCHNG    10
#define TAG_DATA_BUFFEXCHNG       11
#define TAG_NUMDATA_PARTIOEXCHNG  12
#define TAG_DATA_PARTIOEXCHNG     13
#define TAG_NUMDATA1_FOFEXCHNG    14
#define TAG_DATA1_FOFEXCHNG       15
#define TAG_NUMDATA2_FOFEXCHNG    16
#define TAG_DATA2_FOFEXCHNG       17

#define IO_CHUNK_SIZE_MB 100.0
#define EXTRA_PART_ALLOC_FACTOR 1.20

typedef struct {
  float OmegaM;
  float h;
  float PartMass;
  long TotNumParts;
  float ScaleFactorSnapshot;
  float BoxLengthSnapshot;
  float SofteningLength;
  char SimulationType[MAX_FILENAME];
  char HaloOutputFile[MAX_FILENAME];
  char PeakOutputFile[MAX_FILENAME];
  char PeakInputFile[MAX_FILENAME];
  char SnapshotFileList[MAX_FILENAME]; /* file(s) to read particle data from */
  char SnapshotHeaderList[MAX_FILENAME]; /* file(s) to read particle data from */
  int NumTasksIOInParallel; /* the number of tasks which are allowed to read from disk at the same time */
  float PeakApertureSizeFactor; /* size in softening lengths of fixed aperture to find peaks within */
  float MinimumPeakOverdensWRTMeanDens; /* minimum overdensity wrt to mean dens of snapshot for finding intitial peaks */
  float MinNumPartsHalo; /* min number of parts in a halo */
  float DomainBuffSize; /* size of buffer between domains for halo finding */
  char HaloTruncRadType[MAX_FILENAME]; /* sets type of truncation rule for halos */
  float HaloOverdens;  /* for use with overdensity criterion halo defs - compile options HALOMEANOVERDENS and HALOCRITOVERDENS */
  char HaloExclType[MAX_FILENAME]; /* sets type of exclusion rule for halos */
} HaloFinderData;

typedef struct {
  float pos[3];
  float vel[3];
  long index;
  float mass;
  float radius;
  float vcirc;
} Halo;

typedef struct {
  int length;
  int link;
  int head;
  int tail;
} FoFData;

//from gadget-2 peano.c functions 
typedef int peanokey;

#define BITS_PER_DIMENSION 10 /*!< Bits per dimension available for Peano-Hilbert order.                                                                       
				Note: If peanokey is defined as type int, the allowed maximum is 10.                                                                                             
				If 64-bit integers are used, the maximum is 21 */
//#define  PEANOCELLS (((peanokey)1)<<(3*BITS_PER_DIMENSION))  /*!< The number of different Peano-Hilbert cells */

/* extern defs of global vars in globalvars.c */
extern HaloFinderData haloFinderData;                /* global struct with all vars from config file */
extern int ThisTask;                                 /* this task's rank in MPI_COMM_WORLD */
extern int NTasks;                                   /* number of tasks in MPI_COMM_WORLD */
extern int NumDomains[3];                            /* # of splits in each direction */
extern int domainCoords[3];                          /* index # of domain in each direction */
extern MPI_Comm cartComm;                            /* MPI_Comm for cartesian topology */

extern long NumPartsInDomain;                        /* number of parts in ThisTask's domain */
extern long NumParts;                                /* number of particles for ThisTask */
extern float *PartsPx;                               /* particle x position */
extern float *PartsPy;                               /* particle y position */
extern float *PartsPz;                               /* particle z position */
extern float *PartsVx;                               /* particle x velocity */
extern float *PartsVy;                               /* particle y velocity */
extern float *PartsVz;                               /* particle z velocity */
extern short int *PartsInDomainMark;                 /* marks if a particle is in the domain or a buffer particle */
extern kdTreeData *PartskdTreeData;                  /* kdTree data struct for particles */

extern long NumHalos;                                /* number of halos for ThisTask */
extern Halo *Halos;                                  /* halos for ThisTask */
extern MPI_Datatype MPI_Halo;                        /* MPI halo data type */

/* in readART.c */
void reorder(char *buffer, size_t size);
int read_ART_Page(long PageNum, long *NumBuffParts, float **buffPx, float **buffPy, float **buffPz, float **buffVx, float **buffVy, float **buffVz);
long read_ART_NumPages(void);

/* in readGADGET.c */
int read_GADGET_Page(long PageNum, long *NumBuffParts, float **buffPx, float **buffPy, float **buffPz, float **buffVx, float **buffVy, float **buffVz);
long read_GADGET_NumPages(void);

/* in readGADGETsinglefile.c */
int read_GADGETsinglefile_Page(long PageNum, long *NumBuffParts, float **buffPx, float **buffPy, float **buffPz, float **buffVx, float **buffVy, float **buffVz);
long read_GADGETsinglefile_NumPages(void);

/* in config.c */
void read_config_file(char fname[]);

/* in partio.c */
void readpartsfromsnapshot(void);
void distributepartsfromsnapshot(float *xx[6], int *sendcounts, int *displs);
void exchangebufferregions(void);

/* in peakhaloasciiio.c */
long fnumlines(FILE *fp);
void writehalos(void);
void readhalos(void);
void writepeaks(void);
void readpeaks(void);

/* in peakhalofitsio.c */
#ifdef HAVEFITSIO
int writepeakshalos_name(char *name);
long read_totnumpeakshalos_fits(char name[]);
long read_chunksizepeakshalos_fits(char name[]);
void read_peakshaloschunk_fits(char name[], Halo *chunkHalos, long BaseHaloChunkSize, long HaloChunkSize, long chunkInd);
#endif

/* in halodefs.c */
float halo_truncation_radius(double m);
float halo_exlcusion_radius(double m1, double m2);

/* in fortran_io_utils.c */
size_t fwrite_fortran(const void *ptr, size_t size, size_t nitems, FILE *stream);
size_t fread_fortran(void *ptr, size_t size, size_t nitems, FILE *stream);

/* halofinder_utils.c */
void destroy_parts_kdtree(void);
void build_parts_kdtree(void);
void reorder_parts_tree(void);
void destroyhalos(void);
void reallocparts(long NumExtra);
void destroyparts(void);
void wrapdist(float *dx, float rcube, float rhalf);
void wrappos(float *dx, float rbox);
peanokey peano_hilbert_key(int x, int y, int z, int bits);
void sort_halos_vcirc(void);
void sort_halos_mass(void);

/* in domain.c */
void getDomainOrigin_shiftpartshalos(float domainOrigin[3]);
void point_global2domain(float *px, float *py, float *pz, float domainOrigin[3]);
void point_domain2global(float *px, float *py, float *pz, float domainOrigin[3]);
void setupdomaindecomp(void);
void destroydomaindecomp(void);
void shiftpartshalosglobal2domain(void);
void shiftpartshalosdomain2global(void);
void shifthalosglobal2domain(void);
void shifthalosdomain2global(void);
void scatterpeakshalos(void);
void gatherpeakshalos(void);

/* in halofinder.c */
void computehaloradii(void);

/* in halopercolation_serial.c */
void halopercolation_serial(void);

/* in halopercolation_parallel.c */
void halopercolation_parallel(void);
void percolate_fofhalogroups(int periodic, FoFData *fofd);
void join_fofhalogroups(kdTreeData **td, float **px, float **py, float **pz, int periodic, float globalMaxHaloMass,
                        float **nbrsRad2, int **nbrsInd, int *maxNumNbrs, FoFData **fofd);
FoFData *build_fofhalogroups(kdTreeData *td, float *px, float *py, float *pz, int periodic,
                             float globalMaxHaloMass, float **nbrsRad2, int **nbrsInd, int *maxNumNbrs);

/* in peakfinder.c */
void peakfinder(void);
void peakpercolation(void);


#ifdef TEST_CODE
/* in test_code.c */
void test_fofdata(FoFData *fofd);
void test_fofhalogroupfinder(void);
void test_sort_halos(void);
void test_kdtree_knnbrsfind(void);
void test_peakhaloio(void);
void test_exchangebufferregions(void);
void test_shiftpartshalos(void);
void test_scattergatherpeakshalos(void);
void test_kdtree_nbrsfind(void);
#endif

#endif /* _HALOFINDER_ */
